"""
file            pylinpro.py
author      ernesto.adorio@gmail.com
                UPEPP UP at Clark Field
                Pampanga, the PHilippines 
                zip code: 2009
source      translation  from BASIC to C/C++ to Python
                of G.E. linpro program.
revisions   2009.01.25  first version
                2009.01.26  debugged and updated.
"""
 
from matlib import matzero, matdim, matprint, matsetblock
 
def pylinpro (AA,M,N,L,  E,G, artifcode,opticode,Basis, disp_code, ZTOL = 1.0e-6):
  """
  Parameters
  AA           AA is a dynamically allocated matrix
  and should be dimensioned (M+2) * JLIM
  where
  JLIM = N + G + L + 1 if artifcode = 0
  and
  JLIM = N + G + M + 1 if artifcode = 1.
 
  B coefficients are stored in AA[ ][0]
  Z coefficients are stored in AA[M][]
  Phase 1 coefficients are stored in AA[M+1][]
 
  M             Number of constraint equations
  N             Number of variables
  L             Number of less thans
  E             Number of equal tos
  G             Number of greater thans
 
  artifcode
  0 - do not process artificial variables,
  1 - otherwise
 
  opticode
  0 - minimization
  1 - maximization
 
  Basis[]  Output current Basis column.
  Variables are numbered 1 to L+E+G.
 
  LPcode   Return LPcode code of LP routine:
  0 - Successful
  1 - Array too small
  2 - Inconsistent equations
  3 - B coefficients must be positive !
  4 - Infeasible, artificial variable in solution
  5 - No solution, unbounded solution
  6 - Excessive number of iterations
  7 - Failed to allocate space for VC
 
  disp_code
  0   no display
  1   for display of results only
  2   for display of Basis and objective function
  3   for display of tableau at each iteration
 
  Description
  This is a C++ translation of the G.E.(General Electric)
  LINPRO program.
  """
 
  def disp_Basis():
      print "Iteration No. %d " % LPIter
      print "Basis            Value "
      # print(" Variable ")
      for sub_i in range(M): # (sub_i = 0 sub_i < M sub_i++):
          print " %4d  %15.5lf" % (Basis[sub_i], AA[sub_i][0])
 
      print "The value of the Objective function is %11.5lf" % (-AA[M][0] if (opticode == 0) else AA[M][0])#? -AA[M][0] : AA[M][0])
 
  def disp_Tableau():
      print "Tableau at iteration : %-d " % LPIter
      print "Basis      Value        ",
      for sub_j in range(1, JLIM): # (sub_j = 1 sub_j < JLIM sub_j++):
          print "%-8d" % VC[sub_j],
 
      print ""
 
      for sub_i in range(M+2): #(sub_i = 0 sub_i < M + 2 sub_i++):
          if (sub_i == M):
              print "%4c" % ' ',
              for sub_j in range((JLIM+1)*8):  # (sub_j = 0 sub_j < (JLIM + 1) * 8 sub_j++):
                  print "%c" % '=',
              print ""
 
          if (sub_i < M):
              print "%8d" % Basis[sub_i],
 
          else:
              print "%8c" % ' ',
 
 
          print "%12.3lf" % AA[sub_i][0],
          for sub_j in range(1, JLIM):   #(sub_j = 1 sub_j < JLIM sub_j++):
              print "%8.3lf" % AA[sub_i][sub_j],
          print 
 
 
  MAXDOUBLE = 1.0e100
  print "M = %d" % M
  print "N = %d" % N
  print "L = %d" % L
  print "E = %d" % E
  print "G = %d" % G
 
  K = N + G + M + 1
  if (artifcode):
      JLIM = K
  else :
      JLIM = N + G + L + 1
 
 
 
  # Check for  size of AA, and equation for L, E and G 
  nrows, ncols  = matdim(AA)
 
  print "nrows = %d" % nrows
  print "ncols = %d" % ncols
  print "JLIM = %d"  % JLIM
  print "M = %d" % M
 
  if ((nrows < M + 2) or ((ncols < JLIM))):
      LPcode = 1		#  Array too small  
      return LPcode
 
  elif (L + E + G != M):
      LPcode = 2		#  Inconsistent equations. 
      return LPcode
 
 
  # Check for B 
  for I in range(M+1): #(I = 0 I <= M I++):
      if (AA[I][0] < 0.0):
	  LPcode = 3
	  return LPcode
 
 
 
  # Initialize the Basis 
  for I in range(M): # (I = 0 I < M I++):
      Basis[I] = I + 1 + N + G
 
  for I in range(M, M+2): #(I = M I <= M + 1 I++):
      Basis[I] = 0
 
 
 
  # Initialize the Variable Header  
  VC = range(JLIM) #(int *) (malloc (sizeof (int) * JLIM))
  if (VC is None):
      LPcode = 7
      return LPcode
 
 
  for J in range(1, JLIM):# (J = 1 J < JLIM J++):
      VC[J] = J
 
 
  # Initialize to zero the surplus and slack regions
  for I in range(M+2): #(I = 0 I <= M + 1 I++):
      for J in range(N+1, JLIM): #(J = N + 1 J < JLIM J++):
	  AA[I][J] = 0.0
 
 
 
  # The  surplus part
  J = N + 1
  for I in range(L + E, M): #(I = L + E I < M I++):
      AA[I][J] = -1.0
      J+= 1
 
 
  # The slack part of the array (include artificial
  # part if JLIM == K)
  J = N + G + 1
  for I in range(M): #(I = 0 I < M I++):
      AA[I][J] = 1.0
      J+= 1
      if (J >= JLIM):
	break
 
 
  # Negate the objective coefficients for
  # maximization problems
  if (opticode == 1):
      for J in range(1, JLIM):  #(J = 1 J < JLIM J++):
	  AA[M][J] = -AA[M][J]
 
 
 
  # Compute the last row for phase 1 
  do_phase1 = False if (G + E == 0) else True # 
  if (do_phase1):
      for J in range(JLIM): #(J = 0 J < JLIM J++):
	  RSum = 0.0
	  for I in range(L, M): #(I = L I < M I++):
	      RSum = RSum + AA[I][J]
 
	  AA[M + 1][J] = -RSum
 
 
  else :
      for J in range(JLIM): #(J = 0 J < JLIM J++):
	  AA[M + 1][J] = 0.0
 
 
  Test_Row = M+1 if (do_phase1) else M # ? M + 1 : M
  maxLPIter = 2 * K
  for LPIter in range(1, maxLPIter+1): #(LPIter = 1 LPIter <= maxLPIter LPIter++):
 
      # Step 1:  find pivot column 
      while (1):
	  if (disp_code == 2):
	    disp_Basis()
	  if (disp_code == 3):
	    disp_Tableau()
 
	  # find the most negative number in the test row 
	  tempR = 0.0
	  Piv_Col = 0
	  for J in range(JLIM -1, 0, -1): #(J = JLIM - 1 J >= 1 J--):
	      R = AA[Test_Row][J]
	      if (R < tempR):
		  if (do_phase1 or (AA[Test_Row + 1][J] == 0)):
		      tempR = R
		      Piv_Col = J
 
 
 
 
	  if ((Piv_Col == 0) and do_phase1):
	      Test_Row-= 1
	      do_phase1 = False
	      if (disp_code == 3):
		  print " End of Phase 1 "
 
 
	  else:
	      break
 
 
 
      if (Piv_Col == 0):
 
	  for I in range(M): #(I = 0 I < M I++):
	      if ((Basis[I] >= JLIM) and (abs (AA[I][0]) > ZTOL)):
		  if (disp_code >= 1):
		      print "Artificial variable in solution."
 
		  LPcode = 4
 
		  return LPcode
 
 
 
	  if (disp_code >= 1):
	      disp_Basis()
	      print "Succesful termination "
 
 
	  LPcode = 0		# successful termination 
 
	  return LPcode
 
 
      # Step 2: Find pivot row 
      R = MAXDOUBLE
      Piv_Row = -1
      Tie_Row = K + 2
 
      for I in range(M-1, -1, -1): # (I = M - 1 I >= 0 I--):
	  tempR = AA[I][Piv_Col]
	  if (tempR > ZTOL) :
	      Ratio = abs (AA[I][0] / tempR)
	      if (Ratio < R):
		  R = Ratio
		  Piv_Row = I
 
	      elif (Ratio == R): # break row ties 
		  if (Basis[Tie_Row] > Basis[I]):
		      Piv_Row = I
 
 
	      Tie_Row = Basis[I]
 
 
 
      if (Piv_Row == -1):
	  if (disp_code >= 1):
	      print "Unbounded Solution"
 
	  LPcode = 5
 
	  return LPcode
 
 
      if (disp_code >= 2):
	  print "Pivot Row: %d, Pivot Column: %d " % (Piv_Row, Piv_Col)
 
 
      # Step 3: Reduce the array 
      invPivot = 1.0 / AA[Piv_Row][Piv_Col]
      for I in range(Test_Row + 1): #(I = 0 I <= Test_Row I++):
	  R = AA[I][Piv_Col] * invPivot
	  if (I != Piv_Row):
	      for J in range( JLIM): # (J = 0 J < JLIM J++):
		  if (J != Piv_Col):
		      AA[I][J] = AA[I][J] - AA[Piv_Row][J] * R
		      if (abs (AA[I][J]) <= ZTOL):
			  AA[I][J] = 0.0
 
 
 
	      AA[I][Piv_Col] = 0.0
 
 
      for J in range(JLIM):   # (J = 0 J < JLIM J++):
	  AA[Piv_Row][J] = AA[Piv_Row][J] * invPivot
 
      AA[Piv_Row][Piv_Col] = 1.0
 
      # Step 4: replace old Basis variable 
      VC[Piv_Col] = Basis[Piv_Row]
      Basis[Piv_Row] = Piv_Col
 
 
  if (disp_code >= 1):
      print "Maximum iterations exceeded."
 
  LPcode = 6
 
  return LPcode
 
if __name__ == "__main__":
  M = 2
  N = 2
  L = 1
  E = 0
  G = 1
 
  artifcode = False
  opcode = 0
  if (artifcode):
      JLIM = N + G + M + 1
 
  else:
      JLIM = N + G + L + 1
  displaycode = 0
 
  Basis = range(JLIM)
  AA = matzero(M + 2, JLIM)
 
  #matsetblock (AA, 0, 0, M, N,
  #         # problem 1-7, pp.7-8 of Bronson[1982] 
  #        [100000.0, 1.0, 1.0, 0.0, 0.0,
  #         20000.0, 0.0, 0.0, 1.0, 1.0,
  #         40000.0, 1.0, 0.0, 1.0, 0.0,
  #         60000.0, 0.0, 1.0, 0.0, 1.0,
  #         0.0, 1.0, -10.0, 0.0, 0.0,
  #         0.0, 0.0, 0.0, 6.0, -5.0,
  #         0.0, 2.0, -8.0, 0.0, 0.0,
  #         0.0, 0.0, 0.0, 2.0, -8.0,
  #         50000.0, 1.0, 1.0, 0.0, 0.0,
  #         5000.0, 0.0, 0.0, 1.0, 1.0, 
  #         0.0, 4.0, -3.0, 6.0, -1.0])
  matsetblock (AA, 0, 0, M, N,
	       # problem 1-7, pp.7-8 of Bronson[1982] 
	      [10.0, 2.0, 1.0,
	       11.0, 1.0, 2.0,
           0.0, 1.0, 1.0])
  #print "Input matrix AA"
  #matprint(AA)
 
  LPcode = pylinpro (AA, M, N,		# 
	   L, E, G,		# 
	   artifcode,		# process artificial vars?
	   opcode,		# maximization problem
	   Basis, displaycode)
 
  print " LP routine exit code:%d " % LPcode

